Antelope Project Analysis

Behavioral Data

Demographics

First, let’s load the data and transform it into a Wide format table

behav <- read_tsv("data/behav/behav_data.txt",
                  col_types=cols())

behav$material[behav$material=="Swahili"] <- "Vocabulary"
behav$material[behav$material=="US Maps"] <- "Maps"

behav <- behav %>% rename(Gender = gender, Age = age)

wbehav <- behav %>% pivot_wider(values_from = alpha,
                                id_cols = c(Subject, Gender, Age),
                                names_from = material)

At this point, we can look at the participant demographocs

ggplot(data=wbehav, aes(x=Age, col=Gender)) +
  geom_histogram(aes(col="white", fill=Gender),
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 1) +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  ggtitle("Age Distribution by Gender") + 
  theme_pander()

## Warning: Removed 2 rows containing non-finite values (stat_bin).

Rate of Forgetting (Alpha parameter) for the two tasks

We can now examine the distribution and correlation of the rate of forgetting for the verbal (Vocabulary) and Maps:

mu <- behav %>% 
  group_by(material) %>% 
  summarize(alpha=mean(alpha))

## `summarise()` ungrouping output (override with `.groups` argument)

p1 <- ggplot(behav, aes(x=alpha, fill=material)) +
  geom_histogram(col="white", binwidth = 0.025, 
                 alpha=0.5, position="identity") +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  geom_vline(data=mu, aes(xintercept=alpha, color=material),
             linetype="dashed") +
  xlab(expression(paste("Estimated value of ", alpha))) +
  theme_pander() +
  ylab("Count") +
  ggtitle("(A) Distribution By Materials") +
  theme(legend.position = c(0.2, 0.8)) +
  theme(plot.title = element_text(hjust = 0.5))

p2 <- ggplot(wbehav, aes(x=Vocabulary, y=Maps)) +
  geom_point(size=4, alpha=0.5, col="black") + #col=K[7]) +
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() + 
  scale_y_continuous() + 
  ggtitle("(B) Correlation Across Materials") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab(expression(paste(alpha, " Maps"))) +
  geom_rug(col="black", lwd=1, alpha=.5) +
  annotate("segment", x=0.1, y=0.1, xend=0.5, 
           yend=0.5, col="grey", lwd=1, lty=2) +
  theme(plot.title = element_text(hjust = 0.5))
  
grid.arrange(p1, p2, ncol=2)

p3 <- ggplot(wbehav, aes(x=Vocabulary, y=Maps)) +
  geom_point(size=4, alpha=0.5, col="black") + #col=K[7]) +
  geom_smooth(method = "lm", formula = y ~ x, col="red", 
              fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() + 
  scale_y_continuous() + 
  coord_fixed() +
  ggtitle("Correlation of Alpha Parameters\nAcross Materials") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab(expression(paste(alpha, " Maps"))) +
  annotate("segment", x=0.1, y=0.1, xend=0.5, 
           yend=0.5, col="grey", lwd=1, lty=2) +
  annotate("text", label = paste("R =", round(cor(wbehav$Vocabulary,
                                                  wbehav$Maps),
                                              2)), 
                                 x=0.2, y=0.45, col="red") +
  
  theme(plot.title = element_text(hjust = 0.5))

ggMarginal(p3, col="white", fill="black", alpha=0.5,
           type="hist", bins=13)

And here is a summary of the data:

summary(wbehav) %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Subject	Gender	Age	Vocabulary	Maps
X	Length:50	Length:50	Min. :18.00	Min. :0.1219	Min. :0.2120
X.1	Class :character	Class :character	1st Qu.:18.00	1st Qu.:0.2647	1st Qu.:0.3039
X.2	Mode :character	Mode :character	Median :19.00	Median :0.2917	Median :0.3447
X.3	NA	NA	Mean :18.92	Mean :0.2825	Mean :0.3355
X.4	NA	NA	3rd Qu.:19.00	3rd Qu.:0.3203	3rd Qu.:0.3680
X.5	NA	NA	Max. :22.00	Max. :0.3852	Max. :0.4701
X.6	NA	NA	NA’s :2	NA	NA

Eyes-Closed EEG Data

Now we look at the neural correlates of the rates of forgetting in EEG data. We begin with eyes closed EEG data, which is the most common type of resting-state EEG recording.

(Optional) Preprocess the Data

Set the process_raw variable to TRUE to reprocess the raw EEG data. That might take a very long time.

process_raw = F
if (process_raw) {
  setwd("./data/eyes_closed/")
  for (sub in dir()[grep("A[1-9]", dir())]) {
    setwd(sub)
    analyze.logfile(sub, "closed")
    setwd("..")
  }
  setwd("../..")
}

Load the Eyes-Closed Spectra and Summary Files

The preprocessing step produces a number of text files. Here, we are interested in the summary file, which contains many useful statistics, including the Individual Alpha Frequency (i.e., the Alpha Peak) for each channel, the spectra file, which contains the estimated log power for each channel in increments of 0.5 Hz.

ec_spectra <- NULL
for (sub in dir("data/eyes_closed/")[grep("A[1-9]",
                                        dir("data/eyes_closed/"))]) {
  table <- read_tsv(paste("data/eyes_closed/", sub, "/", 
                          sub, "_closed_spectra.txt", sep=""),
                    col_types = cols())
  if (is.null(ec_spectra)) {
    ec_spectra <- table
  } else {
    ec_spectra <- ec_spectra %>% bind_rows(table)
  }
}

ec_summary <- NULL
for (sub in dir("data/eyes_closed/")[grep("A[1-9]", 
                                          dir("data/eyes_closed/"))]) {
  table <- read_tsv(paste("data/eyes_closed/", sub, "/", 
                          sub, "_closed_summary.txt", sep=""),
                    col_types = cols())
  if (is.null(ec_summary)) {
    ec_summary <- table
  } else {
    ec_summary <- ec_summary %>% bind_rows(table)
  }
}

Individual Alpha frquency

The IAF is perhaps the most important subject level characteristic of the resting state spectrum. For this reason, we want to get some basic statistics. First, we calculate the subject-level IAF by computing the mode IAF across all channels for a given participant.

iaf_cols <- paste(c("AF3", "AF4", "F3", "F4", 
                    "F7", "F8", "FC5", "FC6", 
                    "T7", "T8", "P7", "P8", 
                    "O1", "O2"), 
                  "IAF", 
                  sep="_")

ec_iafs <- ec_summary %>%
  select(c("Subject", iaf_cols)) %>%
  filter(Subject %in% wbehav$Subject) %>%
  pivot_longer(cols = iaf_cols, names_sep="_", 
               names_to = c("Channel", "Discard"), 
               values_to = "IAF_EC") %>%
  select(Subject, Channel, IAF_EC) %>%
  group_by(Subject) %>%
  summarize(IAF_EC = Mode(IAF_EC))

## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(iaf_cols)` instead of `iaf_cols` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.

## `summarise()` ungrouping output (override with `.groups` argument)

The mean IAF thus calculated is 10.3, which is exactly in the middle of our alpha frequency band.

ggplot(data=ec_iafs, aes(x = IAF_EC)) +
  geom_histogram(aes(col="white"), fill="black", 
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 0.5) +
  ggtitle("IAF Distribution, Eyes Closed") + 
  xlab("IAF (Eyes Closed)") +
  theme_pander()

Average Spectrum

First, let’s examine the spectrograms for each channel to make sure it looks normal and the IAF (individual alpha frequency) is reasonable, i.e., the same for each channel and roughly in the middle of the band definition

l_ec_spectra <- pivot_longer(ec_spectra, cols=names(ec_spectra)[3:130], 
                         names_to="Freq")

l_ec_spectra <- l_ec_spectra %>% 
  rename(Power = value) %>%
  add_column(Recording = "Eyes Closed")

l_ec_spectra$Frequency <- as.numeric(substr(l_ec_spectra$Freq, 
                                            0, str_length(l_ec_spectra$Freq) -2))
l_ec_spectra <- l_ec_spectra %>% 
  add_column(Band="Delta", BandMin=0, BandMax=4) 
l_ec_spectra$Band[l_ec_spectra$Frequency <= 40] <- "Gamma"
l_ec_spectra$Band[l_ec_spectra$Frequency < 30] <- "High Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 18] <- "Upper Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 15] <- "Low Beta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 13] <- "Alpha"
l_ec_spectra$Band[l_ec_spectra$Frequency < 8] <- "Theta"
l_ec_spectra$Band[l_ec_spectra$Frequency < 4] <- "Delta"

l_ec_spectra$BandMin[l_ec_spectra$Frequency <= 40] <- 30
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 30] <- 18
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 18] <- 15
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 15] <- 13
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 13] <- 8
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 8] <- 4
l_ec_spectra$BandMin[l_ec_spectra$Frequency < 4] <- 0

l_ec_spectra$BandMax[l_ec_spectra$Frequency <= 40] <- 40
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 30] <- 30
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 18] <- 18
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 15] <- 15
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 13] <- 13
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 8] <- 8
l_ec_spectra$BandMax[l_ec_spectra$Frequency < 4] <- 4

l_ec_spectra$Band <- factor(l_ec_spectra$Band, 
                        levels = c("Delta", "Theta", "Alpha", 
                                   "Low Beta", "Upper Beta",
                                   "High Beta", "Gamma"))

Now, we remove the three participants for which we have poor quality data

l_ec_spectra <- l_ec_spectra %>%
  filter(Subject %in% behav$Subject)

We can visualize the power spectra to ensure that our data looks normal

gd <- l_ec_spectra %>% 
        group_by(Band) %>% 
        summarise(
          Min = mean(BandMin),
          Max = mean(BandMax),
          Power =mean(Power),
          Frequency = mean(BandMin)
          #Channel=
        )

## `summarise()` ungrouping output (override with `.groups` argument)

ggplot(data=l_ec_spectra, aes(x=Frequency, y=Power, Channel)) +
  geom_rect(data = gd, aes(xmin = Min, xmax = Max, fill = Band), 
            ymin=0, ymax=Inf, colour=NA, alpha=0.5) +
  stat_summary(fun.data=mean_sdl, 
               geom = "ribbon", colour = "white", alpha = 0.4) +
  stat_summary(fun = mean, geom = "line", lwd = 1) +
  facet_wrap(~ Channel, ncol=4) +
  scale_alpha_manual(values = seq(0.1, 0.9, 0.1)) +
  ggtitle("Eyes-Closed Power Spectrum Across Channels") +
  ylab("Log Power") + 
  xlab("Frequency (Hz)") + 
  theme_pander() +
  coord_cartesian(xlim=c(1,40), ylim=c(5, 17)) + #ylim=c(5,17)) +
  scale_fill_brewer(palette = "Set3") +
  theme(plot.title = element_text(hjust = 0.5))

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

Now, let’s create a mean spectrogram with annotations that mark the different frequency bands. It serves as a visual reference that our data is not grossly out of whack.

K = brewer.pal(7, "Set3")

al_ec_spectra <- aggregate(l_ec_spectra[c("Power")],
                       list(Subject = l_ec_spectra$Subject, 
                            Frequency = l_ec_spectra$Frequency, 
                            Band = l_ec_spectra$Band), 
                       mean)

ggplot(data=al_ec_spectra, aes(x=Frequency, y=Power, col=Power)) +
  annotate("rect", xmin = 0, xmax = 4, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[1]) +
  annotate("rect", xmin = 4, xmax = 8, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[2]) +
  annotate("rect", xmin = 8, xmax = 13, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[3]) +      
  annotate("rect", xmin = 13, xmax = 15, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[4]) +
  annotate("rect", xmin = 15, xmax = 18, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill = K[5]) +
  annotate("rect", xmin = 18, xmax = 30, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[6]) +
  annotate("rect", xmin = 30, xmax = 40, ymin = 0, ymax = Inf,
           alpha = 0.5, fill=K[7]) +
  stat_summary(fun.data = mean_sdl, geom = "ribbon", 
               alpha = 0.5, fill = "grey50", col = "white") +
  stat_summary(fun=mean, geom="line", lwd=2) +
  coord_cartesian(xlim=c(0,40), ylim=c(5,18)) +
  xlab("Frequency (Hz)") +
  theme_pander() +
  annotate("text", x=2, y=5, label="Delta", angle=90, hjust=0) +
  annotate("text", x=6, y=5, label="Theta", angle=90, hjust=0) +
  annotate("text", x=10.5, y=5, label="Alpha", angle=90, hjust=0) +
  annotate("text", x=14, y=5, label="Low Beta", angle=90, hjust=0) +
  annotate("text", x=16.5, y=5, label="Upper Beta", angle=90, hjust=0) +
  annotate("text", x=24, y=5, label="High Beta", angle=90, hjust=0) +
  annotate("text", x=35, y=5, label="Gamma", angle=90, hjust=0) 

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

Correlations with Rate of Forgetting

Having ensured that our EEG recordings look normal, reflect known neurophysiological characteristics, and that our frequency bands boundaries are reasonable, we can proceed with correlating the rates of forgetting with resting state QEEG characteristics.

First, we calculate the mean power for every band and channel, remove participants who were not included in the behavioral tests, and correct for multiple comparisons by frequency band:

wbehav <- wbehav %>%
  mutate(GlobalAlpha=(Vocabulary + Maps)/2)

Adata_ec <- l_ec_spectra %>% 
  group_by(Subject, Channel, Band, Recording) %>% 
  summarize(Power=mean(Power))

## `summarise()` regrouping output by 'Subject', 'Channel', 'Band' (override with `.groups` argument)

Fdata_ec <- inner_join(Adata_ec, behav, by="Subject")

Rdata_ec <- Fdata_ec %>%
  group_by(material, Channel, Band, Recording) %>%
  summarise(r = cor(Power, alpha),
            p = cor.test(Power, alpha)$p.value)

## `summarise()` regrouping output by 'material', 'Channel', 'Band' (override with `.groups` argument)

Rdata_ec <- Rdata_ec %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH")) 

Here is the full table of statistics:

write_csv(Rdata_ec, "correlations_channel_band_eyes_closed.csv", col_names = T)

Rdata_ec %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

material	Channel	Band	Recording	r	p	q
Maps	AF3	Delta	Eyes Closed	0.010	0.947	0.971
Maps	AF3	Theta	Eyes Closed	0.072	0.621	0.790
Maps	AF3	Alpha	Eyes Closed	0.113	0.436	0.592
Maps	AF3	Low Beta	Eyes Closed	0.189	0.188	0.316
Maps	AF3	Upper Beta	Eyes Closed	0.314	0.026	0.092
Maps	AF3	High Beta	Eyes Closed	0.219	0.127	0.403
Maps	AF3	Gamma	Eyes Closed	0.101	0.484	0.924
Maps	AF4	Delta	Eyes Closed	0.096	0.507	0.971
Maps	AF4	Theta	Eyes Closed	0.079	0.585	0.790
Maps	AF4	Alpha	Eyes Closed	0.147	0.308	0.582
Maps	AF4	Low Beta	Eyes Closed	0.216	0.132	0.316
Maps	AF4	Upper Beta	Eyes Closed	0.339	0.016	0.092
Maps	AF4	High Beta	Eyes Closed	0.280	0.049	0.342
Maps	AF4	Gamma	Eyes Closed	0.145	0.315	0.924
Maps	F3	Delta	Eyes Closed	-0.005	0.971	0.971
Maps	F3	Theta	Eyes Closed	0.072	0.619	0.790
Maps	F3	Alpha	Eyes Closed	0.140	0.333	0.582
Maps	F3	Low Beta	Eyes Closed	0.183	0.203	0.316
Maps	F3	Upper Beta	Eyes Closed	0.256	0.073	0.120
Maps	F3	High Beta	Eyes Closed	0.184	0.200	0.403
Maps	F3	Gamma	Eyes Closed	0.022	0.881	0.924
Maps	F4	Delta	Eyes Closed	0.239	0.095	0.343
Maps	F4	Theta	Eyes Closed	0.122	0.400	0.790
Maps	F4	Alpha	Eyes Closed	0.156	0.278	0.582
Maps	F4	Low Beta	Eyes Closed	0.200	0.164	0.316
Maps	F4	Upper Beta	Eyes Closed	0.290	0.041	0.116
Maps	F4	High Beta	Eyes Closed	0.215	0.133	0.403
Maps	F4	Gamma	Eyes Closed	0.135	0.349	0.924
Maps	F7	Delta	Eyes Closed	-0.058	0.689	0.971
Maps	F7	Theta	Eyes Closed	0.087	0.550	0.790
Maps	F7	Alpha	Eyes Closed	0.115	0.425	0.592
Maps	F7	Low Beta	Eyes Closed	0.203	0.158	0.316
Maps	F7	Upper Beta	Eyes Closed	0.254	0.075	0.120
Maps	F7	High Beta	Eyes Closed	0.121	0.403	0.564
Maps	F7	Gamma	Eyes Closed	-0.021	0.885	0.924
Maps	F8	Delta	Eyes Closed	0.011	0.939	0.971
Maps	F8	Theta	Eyes Closed	0.147	0.309	0.790
Maps	F8	Alpha	Eyes Closed	0.173	0.230	0.582
Maps	F8	Low Beta	Eyes Closed	0.206	0.152	0.316
Maps	F8	Upper Beta	Eyes Closed	0.319	0.024	0.092
Maps	F8	High Beta	Eyes Closed	0.184	0.201	0.403
Maps	F8	Gamma	Eyes Closed	0.065	0.655	0.924
Maps	FC5	Delta	Eyes Closed	-0.030	0.838	0.971
Maps	FC5	Theta	Eyes Closed	0.025	0.865	0.932
Maps	FC5	Alpha	Eyes Closed	0.078	0.591	0.637
Maps	FC5	Low Beta	Eyes Closed	0.128	0.376	0.526
Maps	FC5	Upper Beta	Eyes Closed	0.153	0.288	0.367
Maps	FC5	High Beta	Eyes Closed	0.071	0.622	0.725
Maps	FC5	Gamma	Eyes Closed	0.014	0.924	0.924
Maps	FC6	Delta	Eyes Closed	0.069	0.633	0.971
Maps	FC6	Theta	Eyes Closed	0.106	0.462	0.790
Maps	FC6	Alpha	Eyes Closed	0.167	0.247	0.582
Maps	FC6	Low Beta	Eyes Closed	0.183	0.203	0.316
Maps	FC6	Upper Beta	Eyes Closed	0.272	0.056	0.120
Maps	FC6	High Beta	Eyes Closed	0.153	0.290	0.507
Maps	FC6	Gamma	Eyes Closed	0.088	0.545	0.924
Maps	O1	Delta	Eyes Closed	0.221	0.123	0.343
Maps	O1	Theta	Eyes Closed	0.009	0.952	0.952
Maps	O1	Alpha	Eyes Closed	-0.096	0.508	0.592
Maps	O1	Low Beta	Eyes Closed	-0.039	0.786	0.847
Maps	O1	Upper Beta	Eyes Closed	0.100	0.491	0.573
Maps	O1	High Beta	Eyes Closed	0.093	0.519	0.660
Maps	O1	Gamma	Eyes Closed	0.139	0.336	0.924
Maps	O2	Delta	Eyes Closed	0.229	0.109	0.343
Maps	O2	Theta	Eyes Closed	0.106	0.463	0.790
Maps	O2	Alpha	Eyes Closed	0.060	0.680	0.680
Maps	O2	Low Beta	Eyes Closed	0.098	0.497	0.632
Maps	O2	Upper Beta	Eyes Closed	0.193	0.180	0.252
Maps	O2	High Beta	Eyes Closed	0.205	0.153	0.403
Maps	O2	Gamma	Eyes Closed	0.181	0.209	0.924
Maps	P7	Delta	Eyes Closed	0.201	0.161	0.375
Maps	P7	Theta	Eyes Closed	-0.037	0.798	0.931
Maps	P7	Alpha	Eyes Closed	-0.148	0.306	0.582
Maps	P7	Low Beta	Eyes Closed	-0.028	0.847	0.847
Maps	P7	Upper Beta	Eyes Closed	0.049	0.735	0.791
Maps	P7	High Beta	Eyes Closed	0.027	0.852	0.852
Maps	P7	Gamma	Eyes Closed	0.028	0.846	0.924
Maps	P8	Delta	Eyes Closed	0.326	0.021	0.290
Maps	P8	Theta	Eyes Closed	0.343	0.015	0.208
Maps	P8	Alpha	Eyes Closed	0.222	0.122	0.582
Maps	P8	Low Beta	Eyes Closed	0.298	0.036	0.316
Maps	P8	Upper Beta	Eyes Closed	0.365	0.009	0.092
Maps	P8	High Beta	Eyes Closed	0.332	0.018	0.258
Maps	P8	Gamma	Eyes Closed	0.288	0.043	0.599
Maps	T7	Delta	Eyes Closed	-0.042	0.770	0.971
Maps	T7	Theta	Eyes Closed	-0.148	0.307	0.790
Maps	T7	Alpha	Eyes Closed	-0.099	0.494	0.592
Maps	T7	Low Beta	Eyes Closed	-0.060	0.681	0.795
Maps	T7	Upper Beta	Eyes Closed	-0.023	0.876	0.876
Maps	T7	High Beta	Eyes Closed	-0.037	0.797	0.852
Maps	T7	Gamma	Eyes Closed	-0.017	0.905	0.924
Maps	T8	Delta	Eyes Closed	0.247	0.084	0.343
Maps	T8	Theta	Eyes Closed	0.251	0.078	0.547
Maps	T8	Alpha	Eyes Closed	0.239	0.095	0.582
Maps	T8	Low Beta	Eyes Closed	0.214	0.135	0.316
Maps	T8	Upper Beta	Eyes Closed	0.252	0.077	0.120
Maps	T8	High Beta	Eyes Closed	0.140	0.332	0.517
Maps	T8	Gamma	Eyes Closed	0.036	0.804	0.924
Vocabulary	AF3	Delta	Eyes Closed	-0.008	0.958	0.958
Vocabulary	AF3	Theta	Eyes Closed	0.023	0.873	0.962
Vocabulary	AF3	Alpha	Eyes Closed	0.086	0.553	0.779
Vocabulary	AF3	Low Beta	Eyes Closed	0.298	0.035	0.145
Vocabulary	AF3	Upper Beta	Eyes Closed	0.225	0.116	0.334
Vocabulary	AF3	High Beta	Eyes Closed	0.020	0.890	0.958
Vocabulary	AF3	Gamma	Eyes Closed	-0.091	0.529	0.617
Vocabulary	AF4	Delta	Eyes Closed	-0.059	0.683	0.939
Vocabulary	AF4	Theta	Eyes Closed	0.036	0.802	0.962
Vocabulary	AF4	Alpha	Eyes Closed	0.105	0.466	0.779
Vocabulary	AF4	Low Beta	Eyes Closed	0.257	0.071	0.157
Vocabulary	AF4	Upper Beta	Eyes Closed	0.223	0.119	0.334
Vocabulary	AF4	High Beta	Eyes Closed	0.026	0.859	0.958
Vocabulary	AF4	Gamma	Eyes Closed	-0.114	0.432	0.605
Vocabulary	F3	Delta	Eyes Closed	-0.201	0.162	0.939
Vocabulary	F3	Theta	Eyes Closed	-0.007	0.962	0.962
Vocabulary	F3	Alpha	Eyes Closed	0.064	0.660	0.840
Vocabulary	F3	Low Beta	Eyes Closed	0.235	0.101	0.157
Vocabulary	F3	Upper Beta	Eyes Closed	0.175	0.223	0.428
Vocabulary	F3	High Beta	Eyes Closed	-0.035	0.808	0.958
Vocabulary	F3	Gamma	Eyes Closed	-0.191	0.185	0.581
Vocabulary	F4	Delta	Eyes Closed	-0.041	0.778	0.939
Vocabulary	F4	Theta	Eyes Closed	0.036	0.802	0.962
Vocabulary	F4	Alpha	Eyes Closed	0.105	0.469	0.779
Vocabulary	F4	Low Beta	Eyes Closed	0.273	0.055	0.154
Vocabulary	F4	Upper Beta	Eyes Closed	0.249	0.082	0.334
Vocabulary	F4	High Beta	Eyes Closed	-0.001	0.997	0.997
Vocabulary	F4	Gamma	Eyes Closed	-0.140	0.332	0.581
Vocabulary	F7	Delta	Eyes Closed	0.015	0.917	0.958
Vocabulary	F7	Theta	Eyes Closed	0.029	0.842	0.962
Vocabulary	F7	Alpha	Eyes Closed	0.099	0.494	0.779
Vocabulary	F7	Low Beta	Eyes Closed	0.290	0.041	0.145
Vocabulary	F7	Upper Beta	Eyes Closed	0.168	0.245	0.428
Vocabulary	F7	High Beta	Eyes Closed	-0.043	0.765	0.958
Vocabulary	F7	Gamma	Eyes Closed	-0.126	0.383	0.595
Vocabulary	F8	Delta	Eyes Closed	-0.126	0.383	0.939
Vocabulary	F8	Theta	Eyes Closed	-0.019	0.897	0.962
Vocabulary	F8	Alpha	Eyes Closed	0.137	0.343	0.779
Vocabulary	F8	Low Beta	Eyes Closed	0.241	0.091	0.157
Vocabulary	F8	Upper Beta	Eyes Closed	0.196	0.173	0.403
Vocabulary	F8	High Beta	Eyes Closed	-0.082	0.571	0.958
Vocabulary	F8	Gamma	Eyes Closed	-0.206	0.151	0.581
Vocabulary	FC5	Delta	Eyes Closed	-0.087	0.546	0.939
Vocabulary	FC5	Theta	Eyes Closed	-0.090	0.532	0.962
Vocabulary	FC5	Alpha	Eyes Closed	-0.007	0.962	0.962
Vocabulary	FC5	Low Beta	Eyes Closed	0.148	0.306	0.357
Vocabulary	FC5	Upper Beta	Eyes Closed	0.046	0.751	0.809
Vocabulary	FC5	High Beta	Eyes Closed	-0.115	0.428	0.958
Vocabulary	FC5	Gamma	Eyes Closed	-0.150	0.297	0.581
Vocabulary	FC6	Delta	Eyes Closed	-0.092	0.524	0.939
Vocabulary	FC6	Theta	Eyes Closed	-0.040	0.781	0.962
Vocabulary	FC6	Alpha	Eyes Closed	0.097	0.501	0.779
Vocabulary	FC6	Low Beta	Eyes Closed	0.203	0.157	0.214
Vocabulary	FC6	Upper Beta	Eyes Closed	0.154	0.286	0.444
Vocabulary	FC6	High Beta	Eyes Closed	-0.084	0.562	0.958
Vocabulary	FC6	Gamma	Eyes Closed	-0.172	0.232	0.581
Vocabulary	O1	Delta	Eyes Closed	-0.044	0.759	0.939
Vocabulary	O1	Theta	Eyes Closed	0.057	0.697	0.962
Vocabulary	O1	Alpha	Eyes Closed	0.041	0.776	0.866
Vocabulary	O1	Low Beta	Eyes Closed	0.198	0.168	0.214
Vocabulary	O1	Upper Beta	Eyes Closed	0.138	0.340	0.476
Vocabulary	O1	High Beta	Eyes Closed	-0.052	0.722	0.958
Vocabulary	O1	Gamma	Eyes Closed	-0.145	0.315	0.581
Vocabulary	O2	Delta	Eyes Closed	-0.036	0.804	0.939
Vocabulary	O2	Theta	Eyes Closed	0.077	0.597	0.962
Vocabulary	O2	Alpha	Eyes Closed	0.122	0.399	0.779
Vocabulary	O2	Low Beta	Eyes Closed	0.340	0.016	0.111
Vocabulary	O2	Upper Beta	Eyes Closed	0.271	0.057	0.334
Vocabulary	O2	High Beta	Eyes Closed	0.071	0.625	0.958
Vocabulary	O2	Gamma	Eyes Closed	-0.092	0.524	0.617
Vocabulary	P7	Delta	Eyes Closed	-0.065	0.656	0.939
Vocabulary	P7	Theta	Eyes Closed	-0.121	0.404	0.962
Vocabulary	P7	Alpha	Eyes Closed	-0.085	0.556	0.779
Vocabulary	P7	Low Beta	Eyes Closed	0.040	0.783	0.843
Vocabulary	P7	Upper Beta	Eyes Closed	-0.020	0.892	0.892
Vocabulary	P7	High Beta	Eyes Closed	-0.119	0.410	0.958
Vocabulary	P7	Gamma	Eyes Closed	-0.184	0.201	0.581
Vocabulary	P8	Delta	Eyes Closed	0.072	0.620	0.939
Vocabulary	P8	Theta	Eyes Closed	0.226	0.115	0.962
Vocabulary	P8	Alpha	Eyes Closed	0.210	0.142	0.779
Vocabulary	P8	Low Beta	Eyes Closed	0.424	0.002	0.030
Vocabulary	P8	Upper Beta	Eyes Closed	0.335	0.017	0.242
Vocabulary	P8	High Beta	Eyes Closed	0.158	0.273	0.958
Vocabulary	P8	Gamma	Eyes Closed	0.013	0.926	0.926
Vocabulary	T7	Delta	Eyes Closed	-0.051	0.726	0.939
Vocabulary	T7	Theta	Eyes Closed	-0.104	0.471	0.962
Vocabulary	T7	Alpha	Eyes Closed	0.036	0.804	0.866
Vocabulary	T7	Low Beta	Eyes Closed	-0.027	0.851	0.851
Vocabulary	T7	Upper Beta	Eyes Closed	-0.059	0.684	0.798
Vocabulary	T7	High Beta	Eyes Closed	-0.058	0.690	0.958
Vocabulary	T7	Gamma	Eyes Closed	-0.064	0.660	0.710
Vocabulary	T8	Delta	Eyes Closed	-0.037	0.798	0.939
Vocabulary	T8	Theta	Eyes Closed	0.041	0.779	0.962
Vocabulary	T8	Alpha	Eyes Closed	0.173	0.230	0.779
Vocabulary	T8	Low Beta	Eyes Closed	0.238	0.096	0.157
Vocabulary	T8	Upper Beta	Eyes Closed	0.097	0.504	0.641
Vocabulary	T8	High Beta	Eyes Closed	-0.096	0.506	0.958
Vocabulary	T8	Gamma	Eyes Closed	-0.225	0.116	0.581

And here is the corresponding distributions of r and p values. The dashed lines correspond to a significant threshold of p < .05 on either a r or a p-value scale.

ggplot(Rdata_ec, aes(x = Band, y = r, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  facet_wrap(~ material) +
  ggtitle("Correlation by Frequency Band") +
  ylab("r value") +
  annotate("segment", x=-Inf, xend=Inf, y=0.28, yend = 0.28, lty=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-0.28, yend = -0.28, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme_pander()

ggplot(Rdata_ec, aes(x = Band, y = p, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", alpha=0.5, geom = "errorbar") +
  ggtitle("p-value, by Band") +
  ylab("p value") +
  facet_wrap(~ material) +
  annotate("segment", x=-Inf, xend=Inf, y=0.05, yend = 0.05, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_y_log10() +
  theme_pander()

From the original list, we can extract only those channels that survive the FDR correction:

survivors <- Rdata_ec %>%
  filter(q < 0.05)
  
survivors %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

material	Channel	Band	Recording	r	p	q
Vocabulary	P8	Low Beta	Eyes Closed	0.424	0.002	0.03

Only one channel P8, in one frequency band (low beta) survives FDR correction. We can plot the correlation between power in P8 and rate of forgetting to see the relationship:

focus <- Fdata_ec %>%
  filter(material == "Vocabulary", 
         Channel %in% survivors$Channel,
         Band == "Low Beta") %>%
  rename(Alpha = alpha)

p <- ggplot(focus, aes(x = Alpha, y = Power)) +
  geom_point(size = 4, alpha = 0.5, col = "black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() +   
  scale_y_continuous() + 
  ggtitle("Rate of Forgetting\nand Eyes-Closed Beta Power") +
  xlab(expression(paste(alpha, " Vocabulary"))) +
  ylab("Low Beta (13-15 Hz) Power Over P8") +
  geom_text(data=survivors, col="red",
            mapping=aes(x=0.25, y=10.5, 
                        label= paste("r =", round(r, 2)))) +
  theme(plot.title = element_text(hjust = 0.5))

ggMarginal(p, col = "white", fill = "black", alpha = 0.5,
           type = "hist", bins = 13)

Eyes-Open EEG Data

(Optional) Preprocess the Data

Set the process_raw variable to TRUE to reprocess the raw EEG data. That might take a very long time.

process_raw = F
if (process_raw) {
  setwd("./data/eyes_open/")
  for (sub in dir()[grep("A[1-9]", dir())]) {
    setwd(sub)
    analyze.logfile(sub, "open")
    setwd("..")
  }
  setwd("../..")
}

Load the Eyes-Open Spectra and Summary Files

As in the case of eyes-closed data, we are going to examine the spectral and IAF characteristics of eyes-open recordings though the summary files produced by preprocessing.

eo_spectra <- NULL
for (sub in dir("data/eyes_open/")[grep("A[1-9]", 
                                        dir("data/eyes_open/"))]) {
  table <- read_tsv(paste("data/eyes_open/", sub, "/", 
                          sub, "_open_spectra.txt", sep=""),
                    col_types = cols())
  if (is.null(eo_spectra)) {
    eo_spectra <- table
  } else {
    eo_spectra <- eo_spectra %>% bind_rows(table)
  }
}

eo_summary <- NULL
for (sub in dir("data/eyes_open/")[grep("A[1-9]", 
                                        dir("data/eyes_open/"))]) {
  table <- read_tsv(paste("data/eyes_open/", sub, "/", 
                          sub, "_open_summary.txt", sep=""),
                    col_types = cols())
  if (is.null(eo_summary)) {
    eo_summary <- table
  } else {
    eo_summary <- eo_summary %>% bind_rows(table)
  }
}

Individual Alpha frquency

Since eyes-open recordings are somewhat less common, we need to run a few sanity checks on them. The first and most obsvious concerns the IAF, and whether the IAFs during eyes-open recordings are similar and correlated to those observed during eyes-closed sessions.

eo_iafs <- eo_summary %>%
  select(c("Subject", iaf_cols)) %>%
  filter(Subject %in% wbehav$Subject) %>%
  pivot_longer(cols = iaf_cols, names_sep="_", 
               names_to = c("Channel", "Discard"), 
               values_to = "IAF_EO") %>%
  select(Subject, Channel, IAF_EO) %>%
  group_by(Subject) %>%
  summarize(IAF_EO = Mode(IAF_EO))

## `summarise()` ungrouping output (override with `.groups` argument)

iafs <- inner_join(eo_iafs, ec_iafs, by="Subject")

The mean IAF thus calculated is 10.33, which is, once more, exactly in the middle of our alpha frequency band.

ggplot(data=eo_iafs, aes(x = IAF_EO)) +
  geom_histogram(aes(col="white"), fill="black", 
                 colour="white", alpha=0.5, 
                 position="identity", binwidth = 0.5) +
  ggtitle("IAF Distribution, Eyes-Open") + 
  xlab("IAF (Eyes Open)") +
  theme_pander()

As it can be seen, there is a significant correlation in the IAFs between the two recordings (p = 10^{-6}). Between recordings, the individual IAF remains the same +/- 0.83.

ggplot(iafs, aes(x=IAF_EO, y=IAF_EC)) +
  geom_count(alpha=0.5, col="black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  xlab("IAF Eyes Open (Hz)") +
  ylab("IAF Eyes Closed (Hz") +
  geom_text(data = summarize(iafs, Mean=mean(IAF_EC)),  # Avoid label overlaps
            x= 9.5, y=11.5, 
            col="red", 
            label=paste("r =", 
                        round(cor(iafs$IAF_EC, iafs$IAF_EO), 2))) +
  ggtitle("Correlation Between Eyes-Closed and Eyes-Open IAF") +
  theme_pander()

Average Spectrum

We can also investigate whether the average spectrum is comparable to that of eyes-closed recordings. To do so, we will calculate the average spectrum across channels, much in the same ways as it was done for eyes-closed data:

l_eo_spectra <- pivot_longer(eo_spectra, cols=names(ec_spectra)[3:130], 
                         names_to="Freq")

l_eo_spectra <- l_eo_spectra %>% 
  rename(Power = value) %>%
  add_column(Recording = "Eyes Open")

l_eo_spectra$Frequency <- as.numeric(substr(l_eo_spectra$Freq, 
                                            0, str_length(l_eo_spectra$Freq) -2))
l_eo_spectra <- l_eo_spectra %>% 
  add_column(Band="Delta", BandMin=0, BandMax=4) 
l_eo_spectra$Band[l_eo_spectra$Frequency <= 40] <- "Gamma"
l_eo_spectra$Band[l_eo_spectra$Frequency < 30] <- "High Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 18] <- "Upper Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 15] <- "Low Beta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 13] <- "Alpha"
l_eo_spectra$Band[l_eo_spectra$Frequency < 8] <- "Theta"
l_eo_spectra$Band[l_eo_spectra$Frequency < 4] <- "Delta"

l_eo_spectra$BandMin[l_eo_spectra$Frequency <= 40] <- 30
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 30] <- 18
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 18] <- 15
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 15] <- 13
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 13] <- 8
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 8] <- 4
l_eo_spectra$BandMin[l_eo_spectra$Frequency < 4] <- 0

l_eo_spectra$BandMax[l_eo_spectra$Frequency <= 40] <- 40
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 30] <- 30
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 18] <- 18
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 15] <- 15
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 13] <- 13
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 8] <- 8
l_eo_spectra$BandMax[l_eo_spectra$Frequency < 4] <- 4

l_eo_spectra$Band <- factor(l_eo_spectra$Band, 
                        levels = c("Delta", "Theta", "Alpha", 
                                   "Low Beta", "Upper Beta",
                                   "High Beta", "Gamma"))

Now, we remove the three participants for which we have poor quality data

l_eo_spectra <- l_eo_spectra %>%
  filter(Subject %in% behav$Subject)

And we can visualize the eyes-open power spectra by channel:

gd <- l_eo_spectra %>% 
        group_by(Band) %>% 
        summarise(
          Min = mean(BandMin),
          Max = mean(BandMax),
          Power =mean(Power),
          Frequency = mean(BandMin)
        )

## `summarise()` ungrouping output (override with `.groups` argument)

ggplot(data=l_eo_spectra, aes(x=Frequency, y=Power, Channel)) +
  geom_rect(data = gd, aes(xmin = Min, xmax = Max, fill = Band), 
            ymin=0, ymax=Inf, colour=NA, alpha=0.5) +
  stat_summary(fun.data=mean_sdl, 
               geom = "ribbon", colour = "white", alpha = 0.4) +
  stat_summary(fun = mean, geom = "line", lwd = 1) +
  facet_wrap(~ Channel, ncol=4) +
  scale_alpha_manual(values = seq(0.1, 0.9, 0.1)) +
  ggtitle("Eyes-Open Power Spectrum Across Channels") +
  ylab("Log Power") + 
  xlab("Frequency (Hz)") + 
  theme_pander() +
  coord_cartesian(xlim=c(1,40), ylim=c(5, 17)) + #ylim=c(5,17)) +
  scale_fill_brewer(palette = "Set3") +
  theme(plot.title = element_text(hjust = 0.5))

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

At a first sight, the eyes-open power spectrum looks remarkably similar, with the obvious (and expected) difference of a less prominent pean in the alpha band. To better compare the two types of recordings, we can visualize the avearage spectra on top of ech other:

l_spectra <- rbind(l_ec_spectra, l_eo_spectra)

al_spectra <- l_spectra %>%
  group_by(Subject, Frequency, Band, Recording) %>%
  summarise(Power = mean(Power))

## `summarise()` regrouping output by 'Subject', 'Frequency', 'Band' (override with `.groups` argument)

al_spectra <- al_spectra %>% ungroup()

ggplot(data=al_spectra, aes(x = Frequency, y = Power, 
                            col = Recording, fill = Recording)) +
  annotate("rect", xmin = 0, xmax = 4, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[1]) +
  annotate("rect", xmin = 4, xmax = 8, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[2]) +
  annotate("rect", xmin = 8, xmax = 13, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[3]) +      
  annotate("rect", xmin = 13, xmax = 15, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[4]) +
  annotate("rect", xmin = 15, xmax = 18, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill = K[5]) +
  annotate("rect", xmin = 18, xmax = 30, ymin = 0, ymax = Inf, 
           alpha = 0.5, fill=K[6]) +
  annotate("rect", xmin = 30, xmax = 40, ymin = 0, ymax = Inf,
           alpha = 0.5, fill=K[7]) +
  stat_summary(fun.data = mean_sdl, geom = "ribbon", 
               alpha = 0.5, col = "white") +
  scale_color_brewer(palette = "Set1") +
  scale_fill_brewer(palette = "Set1") +
  stat_summary(fun = mean, geom = "line", lwd = 2) +
  coord_cartesian(xlim=c(0,40), ylim=c(5,18)) +
  xlab("Frequency (Hz)") +
  ylab("Log Power") +
  theme_pander() +
  annotate("text", x=2, y=18, label="Delta", angle=90, hjust=1) +
  annotate("text", x=6, y=18, label="Theta", angle=90, hjust=1) +
  annotate("text", x=10.5, y=18, label="Alpha", angle=90, hjust=1) +
  annotate("text", x=14, y=18, label="Low Beta", angle=90, hjust=1) +
  annotate("text", x=16.5, y=18, label="Upper Beta", angle=90, hjust=1) +
  annotate("text", x=24, y=18, label="High Beta", angle=90, hjust=1) +
  annotate("text", x=35, y=18, label="Gamma", angle=90, hjust=1) 

## Warning: Computation failed in `stat_summary()`:
## Hmisc package required for this function

Correlations with Rate of Forgetting

Like we did for eyes-closed data, we can now calculate the correlations between eyes-open EEG power spectra and the behavioral rate fo forgetting for verbal and visual materials.

wbehav <- wbehav %>%
  mutate(GlobalAlpha=(Vocabulary + Maps) / 2)

Adata_eo <- l_eo_spectra %>% 
  group_by(Subject, Channel, Band, Recording) %>% 
  summarize(Power=mean(Power))

## `summarise()` regrouping output by 'Subject', 'Channel', 'Band' (override with `.groups` argument)

Fdata_eo <- inner_join(Adata_eo, behav, by = "Subject")

Cdata_eo <- Fdata_eo %>%
  group_by(material, Channel, Band, Recording) %>%
  summarise(r = cor(Power, alpha),
            p = cor.test(Power, alpha)$p.value)

## `summarise()` regrouping output by 'material', 'Channel', 'Band' (override with `.groups` argument)

Rdata_eo <- Cdata_eo %>%
  group_by(material, Band, Recording) %>%
  mutate(q = p.adjust(p, method="BH")) 

We save the data so that it can be visualized as a topomap with Python MNE. Here is the full table of statistics:

write_csv(Rdata_eo, "correlations_channel_band_eyes_open.csv", col_names = T)

Rdata_eo %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

material	Channel	Band	Recording	r	p	q
Maps	AF3	Delta	Eyes Open	0.180	0.211	0.686
Maps	AF3	Theta	Eyes Open	0.122	0.397	0.556
Maps	AF3	Alpha	Eyes Open	0.218	0.129	0.200
Maps	AF3	Low Beta	Eyes Open	0.255	0.074	0.206
Maps	AF3	Upper Beta	Eyes Open	0.383	0.006	0.021
Maps	AF3	High Beta	Eyes Open	0.279	0.050	0.140
Maps	AF3	Gamma	Eyes Open	0.236	0.100	0.447
Maps	AF4	Delta	Eyes Open	0.141	0.330	0.686
Maps	AF4	Theta	Eyes Open	0.136	0.347	0.540
Maps	AF4	Alpha	Eyes Open	0.247	0.083	0.184
Maps	AF4	Low Beta	Eyes Open	0.301	0.034	0.206
Maps	AF4	Upper Beta	Eyes Open	0.383	0.006	0.021
Maps	AF4	High Beta	Eyes Open	0.315	0.026	0.140
Maps	AF4	Gamma	Eyes Open	0.256	0.073	0.447
Maps	F3	Delta	Eyes Open	0.065	0.656	0.707
Maps	F3	Theta	Eyes Open	0.246	0.085	0.441
Maps	F3	Alpha	Eyes Open	0.232	0.105	0.184
Maps	F3	Low Beta	Eyes Open	0.221	0.122	0.244
Maps	F3	Upper Beta	Eyes Open	0.326	0.021	0.037
Maps	F3	High Beta	Eyes Open	0.219	0.126	0.253
Maps	F3	Gamma	Eyes Open	0.126	0.384	0.662
Maps	F4	Delta	Eyes Open	0.103	0.478	0.686
Maps	F4	Theta	Eyes Open	0.236	0.099	0.441
Maps	F4	Alpha	Eyes Open	0.248	0.082	0.184
Maps	F4	Low Beta	Eyes Open	0.265	0.063	0.206
Maps	F4	Upper Beta	Eyes Open	0.338	0.016	0.034
Maps	F4	High Beta	Eyes Open	0.294	0.038	0.140
Maps	F4	Gamma	Eyes Open	0.204	0.154	0.447
Maps	F7	Delta	Eyes Open	0.233	0.104	0.686
Maps	F7	Theta	Eyes Open	0.148	0.304	0.540
Maps	F7	Alpha	Eyes Open	0.247	0.084	0.184
Maps	F7	Low Beta	Eyes Open	0.281	0.048	0.206
Maps	F7	Upper Beta	Eyes Open	0.428	0.002	0.013
Maps	F7	High Beta	Eyes Open	0.279	0.050	0.140
Maps	F7	Gamma	Eyes Open	0.202	0.160	0.447
Maps	F8	Delta	Eyes Open	0.071	0.624	0.707
Maps	F8	Theta	Eyes Open	0.175	0.225	0.526
Maps	F8	Alpha	Eyes Open	0.287	0.043	0.184
Maps	F8	Low Beta	Eyes Open	0.242	0.091	0.212
Maps	F8	Upper Beta	Eyes Open	0.336	0.017	0.034
Maps	F8	High Beta	Eyes Open	0.199	0.166	0.291
Maps	F8	Gamma	Eyes Open	0.088	0.544	0.662
Maps	FC5	Delta	Eyes Open	0.032	0.826	0.826
Maps	FC5	Theta	Eyes Open	-0.028	0.845	0.985
Maps	FC5	Alpha	Eyes Open	0.167	0.247	0.315
Maps	FC5	Low Beta	Eyes Open	0.189	0.188	0.263
Maps	FC5	Upper Beta	Eyes Open	0.229	0.109	0.153
Maps	FC5	High Beta	Eyes Open	0.137	0.343	0.403
Maps	FC5	Gamma	Eyes Open	0.094	0.515	0.662
Maps	FC6	Delta	Eyes Open	-0.115	0.425	0.686
Maps	FC6	Theta	Eyes Open	0.094	0.515	0.655
Maps	FC6	Alpha	Eyes Open	0.190	0.187	0.261
Maps	FC6	Low Beta	Eyes Open	0.132	0.362	0.461
Maps	FC6	Upper Beta	Eyes Open	0.179	0.214	0.249
Maps	FC6	High Beta	Eyes Open	0.136	0.346	0.403
Maps	FC6	Gamma	Eyes Open	0.037	0.798	0.798
Maps	O1	Delta	Eyes Open	0.074	0.608	0.707
Maps	O1	Theta	Eyes Open	0.011	0.938	0.988
Maps	O1	Alpha	Eyes Open	0.085	0.556	0.649
Maps	O1	Low Beta	Eyes Open	0.026	0.857	0.857
Maps	O1	Upper Beta	Eyes Open	0.191	0.183	0.233
Maps	O1	High Beta	Eyes Open	0.145	0.314	0.403
Maps	O1	Gamma	Eyes Open	0.126	0.383	0.662
Maps	O2	Delta	Eyes Open	0.123	0.395	0.686
Maps	O2	Theta	Eyes Open	0.137	0.342	0.540
Maps	O2	Alpha	Eyes Open	0.238	0.096	0.184
Maps	O2	Low Beta	Eyes Open	0.210	0.144	0.244
Maps	O2	Upper Beta	Eyes Open	0.344	0.015	0.034
Maps	O2	High Beta	Eyes Open	0.239	0.094	0.220
Maps	O2	Gamma	Eyes Open	0.176	0.222	0.518
Maps	P7	Delta	Eyes Open	0.128	0.374	0.686
Maps	P7	Theta	Eyes Open	-0.002	0.988	0.988
Maps	P7	Alpha	Eyes Open	0.068	0.640	0.689
Maps	P7	Low Beta	Eyes Open	0.041	0.778	0.838
Maps	P7	Upper Beta	Eyes Open	0.112	0.440	0.474
Maps	P7	High Beta	Eyes Open	0.074	0.608	0.608
Maps	P7	Gamma	Eyes Open	0.062	0.669	0.720
Maps	P8	Delta	Eyes Open	0.117	0.419	0.686
Maps	P8	Theta	Eyes Open	0.275	0.053	0.441
Maps	P8	Alpha	Eyes Open	0.316	0.025	0.184
Maps	P8	Low Beta	Eyes Open	0.320	0.024	0.206
Maps	P8	Upper Beta	Eyes Open	0.436	0.002	0.013
Maps	P8	High Beta	Eyes Open	0.311	0.028	0.140
Maps	P8	Gamma	Eyes Open	0.237	0.098	0.447
Maps	T7	Delta	Eyes Open	-0.100	0.490	0.686
Maps	T7	Theta	Eyes Open	-0.191	0.184	0.515
Maps	T7	Alpha	Eyes Open	-0.033	0.818	0.818
Maps	T7	Low Beta	Eyes Open	-0.092	0.524	0.611
Maps	T7	Upper Beta	Eyes Open	-0.072	0.621	0.621
Maps	T7	High Beta	Eyes Open	-0.137	0.344	0.403
Maps	T7	Gamma	Eyes Open	-0.105	0.468	0.662
Maps	T8	Delta	Eyes Open	0.198	0.168	0.686
Maps	T8	Theta	Eyes Open	0.219	0.126	0.441
Maps	T8	Alpha	Eyes Open	0.270	0.058	0.184
Maps	T8	Low Beta	Eyes Open	0.203	0.157	0.244
Maps	T8	Upper Beta	Eyes Open	0.234	0.102	0.153
Maps	T8	High Beta	Eyes Open	0.081	0.576	0.608
Maps	T8	Gamma	Eyes Open	0.083	0.567	0.662
Vocabulary	AF3	Delta	Eyes Open	0.157	0.277	0.894
Vocabulary	AF3	Theta	Eyes Open	0.042	0.770	0.953
Vocabulary	AF3	Alpha	Eyes Open	0.176	0.220	0.343
Vocabulary	AF3	Low Beta	Eyes Open	0.194	0.177	0.394
Vocabulary	AF3	Upper Beta	Eyes Open	0.218	0.128	0.358
Vocabulary	AF3	High Beta	Eyes Open	0.130	0.370	0.991
Vocabulary	AF3	Gamma	Eyes Open	0.080	0.583	0.720
Vocabulary	AF4	Delta	Eyes Open	0.090	0.536	0.894
Vocabulary	AF4	Theta	Eyes Open	0.062	0.666	0.953
Vocabulary	AF4	Alpha	Eyes Open	0.196	0.173	0.343
Vocabulary	AF4	Low Beta	Eyes Open	0.267	0.061	0.282
Vocabulary	AF4	Upper Beta	Eyes Open	0.253	0.077	0.269
Vocabulary	AF4	High Beta	Eyes Open	0.139	0.334	0.991
Vocabulary	AF4	Gamma	Eyes Open	0.079	0.586	0.720
Vocabulary	F3	Delta	Eyes Open	-0.136	0.346	0.894
Vocabulary	F3	Theta	Eyes Open	0.042	0.772	0.953
Vocabulary	F3	Alpha	Eyes Open	0.165	0.253	0.354
Vocabulary	F3	Low Beta	Eyes Open	0.087	0.549	0.641
Vocabulary	F3	Upper Beta	Eyes Open	0.182	0.206	0.380
Vocabulary	F3	High Beta	Eyes Open	-0.002	0.991	0.991
Vocabulary	F3	Gamma	Eyes Open	-0.096	0.509	0.720
Vocabulary	F4	Delta	Eyes Open	-0.122	0.398	0.894
Vocabulary	F4	Theta	Eyes Open	0.043	0.768	0.953
Vocabulary	F4	Alpha	Eyes Open	0.177	0.218	0.343
Vocabulary	F4	Low Beta	Eyes Open	0.176	0.222	0.394
Vocabulary	F4	Upper Beta	Eyes Open	0.201	0.163	0.379
Vocabulary	F4	High Beta	Eyes Open	0.010	0.946	0.991
Vocabulary	F4	Gamma	Eyes Open	-0.114	0.430	0.720
Vocabulary	F7	Delta	Eyes Open	0.156	0.279	0.894
Vocabulary	F7	Theta	Eyes Open	0.087	0.550	0.953
Vocabulary	F7	Alpha	Eyes Open	0.188	0.192	0.343
Vocabulary	F7	Low Beta	Eyes Open	0.215	0.133	0.394
Vocabulary	F7	Upper Beta	Eyes Open	0.258	0.070	0.269
Vocabulary	F7	High Beta	Eyes Open	0.142	0.325	0.991
Vocabulary	F7	Gamma	Eyes Open	0.072	0.617	0.720
Vocabulary	F8	Delta	Eyes Open	-0.050	0.729	0.894
Vocabulary	F8	Theta	Eyes Open	0.001	0.995	0.995
Vocabulary	F8	Alpha	Eyes Open	0.208	0.148	0.343
Vocabulary	F8	Low Beta	Eyes Open	0.175	0.225	0.394
Vocabulary	F8	Upper Beta	Eyes Open	0.178	0.217	0.380
Vocabulary	F8	High Beta	Eyes Open	0.013	0.928	0.991
Vocabulary	F8	Gamma	Eyes Open	-0.104	0.473	0.720
Vocabulary	FC5	Delta	Eyes Open	-0.073	0.617	0.894
Vocabulary	FC5	Theta	Eyes Open	-0.172	0.233	0.953
Vocabulary	FC5	Alpha	Eyes Open	0.063	0.664	0.721
Vocabulary	FC5	Low Beta	Eyes Open	0.049	0.737	0.737
Vocabulary	FC5	Upper Beta	Eyes Open	0.085	0.559	0.711
Vocabulary	FC5	High Beta	Eyes Open	0.006	0.970	0.991
Vocabulary	FC5	Gamma	Eyes Open	-0.019	0.897	0.897
Vocabulary	FC6	Delta	Eyes Open	-0.166	0.248	0.894
Vocabulary	FC6	Theta	Eyes Open	-0.037	0.799	0.953
Vocabulary	FC6	Alpha	Eyes Open	0.110	0.446	0.568
Vocabulary	FC6	Low Beta	Eyes Open	0.109	0.450	0.610
Vocabulary	FC6	Upper Beta	Eyes Open	0.054	0.710	0.764
Vocabulary	FC6	High Beta	Eyes Open	-0.018	0.902	0.991
Vocabulary	FC6	Gamma	Eyes Open	-0.101	0.485	0.720
Vocabulary	O1	Delta	Eyes Open	-0.096	0.506	0.894
Vocabulary	O1	Theta	Eyes Open	-0.049	0.734	0.953
Vocabulary	O1	Alpha	Eyes Open	0.185	0.198	0.343
Vocabulary	O1	Low Beta	Eyes Open	0.102	0.480	0.610
Vocabulary	O1	Upper Beta	Eyes Open	0.163	0.258	0.401
Vocabulary	O1	High Beta	Eyes Open	-0.033	0.818	0.991
Vocabulary	O1	Gamma	Eyes Open	-0.108	0.457	0.720
Vocabulary	O2	Delta	Eyes Open	-0.043	0.766	0.894
Vocabulary	O2	Theta	Eyes Open	0.026	0.859	0.953
Vocabulary	O2	Alpha	Eyes Open	0.242	0.090	0.343
Vocabulary	O2	Low Beta	Eyes Open	0.280	0.049	0.282
Vocabulary	O2	Upper Beta	Eyes Open	0.293	0.039	0.269
Vocabulary	O2	High Beta	Eyes Open	0.058	0.692	0.991
Vocabulary	O2	Gamma	Eyes Open	-0.081	0.574	0.720
Vocabulary	P7	Delta	Eyes Open	-0.123	0.395	0.894
Vocabulary	P7	Theta	Eyes Open	-0.171	0.234	0.953
Vocabulary	P7	Alpha	Eyes Open	0.062	0.670	0.721
Vocabulary	P7	Low Beta	Eyes Open	-0.054	0.712	0.737
Vocabulary	P7	Upper Beta	Eyes Open	0.001	0.994	0.994
Vocabulary	P7	High Beta	Eyes Open	-0.127	0.378	0.991
Vocabulary	P7	Gamma	Eyes Open	-0.155	0.281	0.720
Vocabulary	P8	Delta	Eyes Open	-0.067	0.646	0.894
Vocabulary	P8	Theta	Eyes Open	0.100	0.491	0.953
Vocabulary	P8	Alpha	Eyes Open	0.241	0.092	0.343
Vocabulary	P8	Low Beta	Eyes Open	0.296	0.037	0.282
Vocabulary	P8	Upper Beta	Eyes Open	0.313	0.027	0.269
Vocabulary	P8	High Beta	Eyes Open	0.082	0.574	0.991
Vocabulary	P8	Gamma	Eyes Open	-0.049	0.736	0.793
Vocabulary	T7	Delta	Eyes Open	0.006	0.968	0.968
Vocabulary	T7	Theta	Eyes Open	-0.202	0.160	0.953
Vocabulary	T7	Alpha	Eyes Open	0.012	0.935	0.935
Vocabulary	T7	Low Beta	Eyes Open	-0.109	0.450	0.610
Vocabulary	T7	Upper Beta	Eyes Open	-0.097	0.502	0.703
Vocabulary	T7	High Beta	Eyes Open	-0.132	0.359	0.991
Vocabulary	T7	Gamma	Eyes Open	-0.108	0.454	0.720
Vocabulary	T8	Delta	Eyes Open	0.006	0.968	0.968
Vocabulary	T8	Theta	Eyes Open	0.021	0.885	0.953
Vocabulary	T8	Alpha	Eyes Open	0.179	0.213	0.343
Vocabulary	T8	Low Beta	Eyes Open	0.186	0.196	0.394
Vocabulary	T8	Upper Beta	Eyes Open	0.057	0.693	0.764
Vocabulary	T8	High Beta	Eyes Open	-0.098	0.497	0.991
Vocabulary	T8	Gamma	Eyes Open	-0.103	0.477	0.720

And here is the corresponding distributions of r and p values. The dashed lines correspond to a significant threshold of p < .05 on either a r or a p-value scale.

ggplot(Rdata_eo, aes(x = Band, y = r, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  facet_wrap(~ material) +
  ggtitle("Correlation, by Band") +
  ylab("r value") +
  annotate("segment", x=-Inf, xend=Inf, y=0.28, yend = 0.28, lty=2) +
  annotate("segment", x=-Inf, xend=Inf, y=-0.28, yend = -0.28, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  theme_pander()

ggplot(Rdata_eo, aes(x = Band, y = p, col = Channel)) +
  geom_point() +
  stat_summary(fun.data = "mean_se", col="black", 
               alpha=0.5, geom = "errorbar") +
  ggtitle("p-value, by Band") +
  ylab("p value") +
  facet_wrap(~ material) +
  annotate("segment", x=-Inf, xend=Inf, y=0.05, yend = 0.05, lty=2) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  scale_y_log10() +
  theme_pander()

Now, we might wonder how much the correlations change from eyes closed to eyes open recordings. The best way to chack is to visualize then as a scatterplot:

Rdata <- bind_rows(Rdata_ec, Rdata_eo)
wrdata <- pivot_wider(Rdata, id_cols = c("material", "Channel", "Band"),
                      names_from = Recording, 
                      values_from = r)

ggplot(wrdata, aes(x = `Eyes Open`, y = `Eyes Closed`, 
                   col = material)) +
  annotate("segment", x = 0, y=-Inf, xend=0, 
           yend=Inf, col="grey", lwd=1, lty=1) +
  annotate("segment", x = -Inf, y=0, xend=Inf, 
           yend=0, col="grey", lwd=1, lty=1) +
  
  geom_point(alpha=0.5, size=3) +
  annotate("segment", x = -0.25, y=-0.25, xend=0.5, 
           yend=0.5, col="black", lwd=1, lty=2) +
  theme(plot.title = element_text(hjust = 0.5)) +
  #geom_smooth(method = "lm", formula = y ~ x, 
  #            fullrange = T, lwd = 1) +
  scale_color_brewer(palette="Dark2") +
  ggtitle("Correlation Coefficients Across Recordings") +
  theme_pander()

Let’s tabulate only those channels that survive the FDR correction:

survivors <- Rdata_eo %>%
  filter(q < 0.05)
  
survivors %>%
  xtable() %>%
  kable(digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

material	Channel	Band	Recording	r	p	q
Maps	AF3	Upper Beta	Eyes Open	0.383	0.006	0.021
Maps	AF4	Upper Beta	Eyes Open	0.383	0.006	0.021
Maps	F3	Upper Beta	Eyes Open	0.326	0.021	0.037
Maps	F4	Upper Beta	Eyes Open	0.338	0.016	0.034
Maps	F7	Upper Beta	Eyes Open	0.428	0.002	0.013
Maps	F8	Upper Beta	Eyes Open	0.336	0.017	0.034
Maps	O2	Upper Beta	Eyes Open	0.344	0.015	0.034
Maps	P8	Upper Beta	Eyes Open	0.436	0.002	0.013

In the case of eyes-open data, we have a much larger set of channels that survive correction for multiple comparison. Although these correlations generally reflect the same trends we have seen the eyes-closed data (i.e., significant positive correlations for prefrontal and right parietal sites), they are now stroger for visual (not verbal) materials and centered at at slightly higher frequency (upper beta instead of low beta).

Here is an overview of the scatterplots and correlations:

focus <- Fdata_eo %>%
  filter(material == "Maps", 
         Channel %in% survivors$Channel,
         Band == "Upper Beta") %>%
  rename(Alpha = alpha)

ggplot(focus, aes(x=Alpha, y=Power)) +
  geom_point(size=4, alpha=0.5, col="black") + 
  geom_smooth(method = "lm", formula = y ~ x, 
              col="red", fill="red", fullrange = T, lwd=2) +
  theme_pander() +
  scale_x_continuous() +   
  scale_y_continuous() + 
  ggtitle("Rate of Forgetting and Eyes-Open Beta Power") +
  xlab(expression(paste(alpha, " Maps"))) +
  ylab("Upper Beta Power (15-18 Hz)") +
  facet_wrap(~ Channel, scales="free_y") +
  geom_text(data=survivors, col="red",
            mapping=aes(x=0.25, y=9.5, 
                        label= paste("r =", round(r, 2)))) +
  theme(panel.spacing = unit(1.5, "lines")) +
  theme(plot.title = element_text(hjust = 0.5))

Predicting Rate of Forgetting with Group Lasso

So far, we have shown that neurophysiological signals exists that are correlated with the rate of forgetting. This, however, does not directly answer wheter we could “read” an individual’s rate of forgetting from their EEG data directy.

Group LASSO

To do so, we need to use a ML approach. Specifically, here we will use group LASSO regresso. Group LASSO is a way to combine regression, validation, and feature selection into a single approach.

LASSO works like a normal regression, with the additional constraint that the \(\beta\) regressors need to minimize the sum of the residual sum of squares (RSS) and the quantity \(- \lambda||\beta||_{1}\), where \(||\beta||_{1}\) is the first-order norm (the sum of absolute values of each regressor weight). For \(\lambda = 0\), LASSO reduces to normal regression but, as the value of \(\lambda\) grows, more and more regressors are “pushed” to zero and removed from the pool. In group LASSO, regressors are combined and selected by groups, which is a natural choice when natural groupings exists, as in the cases of frequency bands.

As a demo, we will pilot this approach to eyes-closed recordings and the combined rates of forgetting.an

Data Preparation

First, we need to transform the data appropriately, converting from the long to the wide format, which is better suited for regression:

library(gglasso)
wfdata <- Fdata_ec %>% 
  pivot_wider(id_cols = c("Subject", "Channel", 
                          "Band", "material", "alpha"), 
              names_from = c("Channel", "Band"), 
              values_from = c("Power"))

wide_data_ec <- wfdata %>%
  pivot_wider(names_from = c("material"), 
              values_from = "alpha")

From the wide format data, we are going to define the data that will enter the regression, the matrix of observations X and the vector of variables Y. The matrix of observations X is 50 x 98 matrix, corresponding to 50 participants x 98 possible power values (14 channels x 7 frequency bands). The Y vector is made of 50 observed rates of forgetting. To make the prediction reliable, we will average the rates of forgetting across materials (Maps + Vocabulary)

channel_names <- c("AF3", "AF4", "F3", "F4", 
                   "F7", "F8", "FC5", "FC6", 
                   "T7", "T8", "P7", "P8", 
                   "O1", "O2")


band_names <- c("Delta", "Theta", "Alpha", "Low Beta", 
                "Upper Beta", "High Beta", "Gamma")

by_band_order <- as.vector(
  outer(channel_names, 
        band_names, 
        function(x, y){
          paste(x, y, sep="_")
          }))

# Predictors
X <- wide_data_ec[,c(-1, -100, -101)] %>%
  #select(ends_with("Beta")) %>%
  as.matrix()

X <- X[,by_band_order]
# We need to reorder X by frequency bands


# Maps
Ymaps <- wide_data_ec[,c(101)] %>%
  as.matrix()

# Vocab
Yvoc <- wide_data_ec[,c(101)] %>%
  as.matrix()

# Average across materials to obrtain stable estimates
Y <- (Ymaps + Yvoc)/2

# Let's save the mean (will be halpful later)
uY <- mean(Y)

# Center Y to remove dominant effect of intercept
Y <- Y - uY

Finally, we need to define our set of grouping variables, defined as a vector of numbers that identify each group. We are going to create two such vectors, one for the freqency bands and one for the channels. In this example we are going to use bands as the grouping variables.

bands <- rep(1:7, each=14)
channels <- sort(rep(1:14, 7))

Fitting the LASSO

Now, with that in place, we can fit the LASSO model. The LASSO model is fit much like a normal regression, with an additional partameter \(\lambda\).

Because we do not know a priori which value of \(lambda\) to pick, the canonical solution is to fit fit different models for different values of \(\lambda\).

fit <- gglasso(x = X, y = Y, 
               group = bands,
               loss = 'ls', 
               nlambda=100,
               lambda.factor = 0.001)
coef.mat=fit$beta

Finding the optimal \(\lambda\) with Cross-Validation

Our goal, of course, is generalizability across participants. Although the fit model is important, what we really care is a different measure, that of cross-validation. For cross-validation, we are going to use the n-fold approach, in which the data is randomly split so that 1/n of it is saved for test. (if n is identical to the number of participants, it becomes the “Leave-One-Out” validation scheme). Here, we used the Leave-One-Out method, which guarantees that the results are the same every time this script is ruin (because there only one possible permutation of the data). So, n = 50.

fit.cv <- cv.gglasso(x = X,
                     y = Y,
                     group = bands,
                     pred.loss = "L1",
                     nlambda = 100,
                     lambda.factor = 0.001,
                     nfolds = 50)

Visualizing Validation Results

Here are the fit and cross-validation plots, with the vertical lines marking the points corresponding to \(\lambda_{min}\) (optimal for cross-validation) and
\(\lambda_{1SE}\) (more conservative model). Because we have used the Leave-One-Out validation scheme, we are can use the the more aggressive \(\lambda_{min}\) estimate, since our sample of the folds was exhaustive.

# Best lambda

lmin <- fit.cv$lambda.min  # Optimal lambda, maximizes generalization
l1se <- fit.cv$lambda.1se  # Conservative model, min + 1SE

plot(fit.cv, 
     pch=20,
     lwd=2)

plot(fit, group=T, 
     col = inferno(length(unique(fit$group)) + 1), 
     lwd = 2)

legend(x="topright", 
       legend=c("Delta", "Theta", "Alpha",
                "Low Beta", "Upper Beta", "High Beta",
                "Gamma"),
       bty="n",
       ncol=2,
       col=inferno(length(unique(fit$group)) + 1),
       lty=1, lwd=2)
abline(v=log(lmin), lty=2)
abline(v=log(l1se), lty=2)

Using the Lambda min value, we can now compare predicted vs. observed values.

prediction <- predict(object = fit,
                    newx = X,
                    s = lmin,
                    type='link')

observed <- data.frame(Subject = wide_data_ec$Subject,
                       RoF = as_vector(Y) + uY,
                       Condition = "Observed") %>%
  rename(RoF = Maps)

predicted <- data.frame(Subject = wide_data_ec$Subject,
                       RoF = as_vector(prediction) + uY,
                       Condition = "Predicted") %>%
  rename(RoF = X1)

comparison <- as_tibble(rbind(observed, predicted))

plt = cbind(Y,
            predict(object = fit,
                    newx = X,
                    s = lmin,
                    type='link'))

# matplot(plt,
#         main="Predicted vs Actual",
#         type='b',
#         lwd=2,
#         pch="o",
#         col=viridis(3),
#         xlab="Subject",
#         ylab="Rate of Forgetting (alpha)")

ggplot(comparison, aes(x=Subject, y=RoF), col=Condition) +
  geom_point(aes(col=Condition)) +
  geom_line(alpha=0.5, lty=2) +
  annotate("text", x= 40, y = 0.20, 
           label=(paste("r =", 
                        round(cor(Y, prediction),
                              2)))) +
  theme_pander() +
  ylab("Rate of Forgetting")+
  ggtitle("Cross-Validation:\nPredicted vs. Observed Rates of Forgetting") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6)) 

wcomparison <- comparison %>%
  pivot_wider(names_from=Condition, values_from = RoF)

p <- ggplot(wcomparison, aes(x = Predicted, y = Observed)) +
  geom_point(alpha=0.5, size=4) +
  geom_abline(slope = 1, intercept = 0,
              lwd=1, col="blue", lty=2) +
  annotate("text", x= 0.40, y = 0.25, 
           label=(paste("r =", 
                        round(cor(Y, prediction),
                              2)))) +
  theme_pander() +
  #coord_fixed() +
  coord_equal(ratio=1) +
  ylab("Rate of Forgetting (Observed)") +
  xlab("Rate of Forgetting (Predicted)") +
  ggtitle("Cross-Validation") 

ggMarginal(p, col="white", fill="black", alpha=0.5,
           type="hist", bins=13)

In terms of overall error, the root mean squared error between predicted and observed Rate of Forgetting is RMSE = 0.024

And now, let’s print the full table of beta regression weights (for inspectio and analysis) at the desiored value of \(\lambda_{min}\). We can also save the weights for future analysis.

get_channel <- function(s) {
  str_split(s, "_", 2)[[1]][1]
}

get_band <- function(s) {
  str_split(s, "_", 2)[[1]][2]
}

vget_channel <- Vectorize(get_channel)
vget_band <- Vectorize(get_band)


coef(fit.cv, 
     s=lmin) -> betas
  
bdf <- data.frame("Predictor" = rownames(betas),
                  "Beta" = as.numeric(betas)) %>%
  mutate(Channel = vget_channel(Predictor),
         Band = vget_band(Predictor)) 

bdf %>%
  filter(Beta != 0) %>%
  xtable() %>%
  kable(digits = 5) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Predictor	Beta	Channel	Band
(Intercept)	0.00955	(Intercept)	NA
AF3_Delta	-0.00659	AF3	Delta
AF4_Delta	0.00294	AF4	Delta
F3_Delta	-0.02081	F3	Delta
F4_Delta	0.03104	F4	Delta
F7_Delta	0.01171	F7	Delta
F8_Delta	-0.01349	F8	Delta
FC5_Delta	-0.01473	FC5	Delta
FC6_Delta	-0.01259	FC6	Delta
T7_Delta	-0.00654	T7	Delta
T8_Delta	0.01771	T8	Delta
P7_Delta	0.00053	P7	Delta
P8_Delta	0.02433	P8	Delta
O1_Delta	0.00076	O1	Delta
O2_Delta	0.01153	O2	Delta
AF3_Theta	-0.00014	AF3	Theta
AF4_Theta	-0.00003	AF4	Theta
F3_Theta	-0.00103	F3	Theta
F4_Theta	-0.00013	F4	Theta
F7_Theta	-0.00101	F7	Theta
F8_Theta	-0.00092	F8	Theta
FC5_Theta	-0.00081	FC5	Theta
FC6_Theta	0.00035	FC6	Theta
T7_Theta	-0.00324	T7	Theta
T8_Theta	0.00128	T8	Theta
P7_Theta	-0.00182	P7	Theta
P8_Theta	0.00146	P8	Theta
O1_Theta	0.00001	O1	Theta
O2_Theta	0.00095	O2	Theta
AF3_Alpha	-0.00175	AF3	Alpha
AF4_Alpha	0.00352	AF4	Alpha
F3_Alpha	0.00034	F3	Alpha
F4_Alpha	0.00289	F4	Alpha
F7_Alpha	0.00115	F7	Alpha
F8_Alpha	0.00520	F8	Alpha
FC5_Alpha	0.00101	FC5	Alpha
FC6_Alpha	0.01122	FC6	Alpha
T7_Alpha	-0.00497	T7	Alpha
T8_Alpha	0.01223	T8	Alpha
P7_Alpha	-0.01927	P7	Alpha
P8_Alpha	-0.00379	P8	Alpha
O1_Alpha	-0.00753	O1	Alpha
O2_Alpha	0.00392	O2	Alpha
AF3_Upper Beta	0.00068	AF3	Upper Beta
AF4_Upper Beta	0.00089	AF4	Upper Beta
F3_Upper Beta	0.00009	F3	Upper Beta
F4_Upper Beta	0.00030	F4	Upper Beta
F7_Upper Beta	0.00042	F7	Upper Beta
F8_Upper Beta	0.00064	F8	Upper Beta
FC5_Upper Beta	0.00014	FC5	Upper Beta
FC6_Upper Beta	0.00069	FC6	Upper Beta
T7_Upper Beta	0.00046	T7	Upper Beta
T8_Upper Beta	0.00027	T8	Upper Beta
P7_Upper Beta	0.00023	P7	Upper Beta
P8_Upper Beta	0.00002	P8	Upper Beta
O1_Upper Beta	0.00058	O1	Upper Beta
O2_Upper Beta	0.00066	O2	Upper Beta
AF3_Gamma	0.00144	AF3	Gamma
AF4_Gamma	0.00539	AF4	Gamma
F3_Gamma	-0.01068	F3	Gamma
F4_Gamma	-0.00540	F4	Gamma
F7_Gamma	-0.01039	F7	Gamma
F8_Gamma	-0.00435	F8	Gamma
FC5_Gamma	-0.00148	FC5	Gamma
FC6_Gamma	-0.00273	FC6	Gamma
T7_Gamma	0.00736	T7	Gamma
T8_Gamma	-0.01448	T8	Gamma
P7_Gamma	-0.00609	P7	Gamma
P8_Gamma	-0.00418	P8	Gamma
O1_Gamma	0.00834	O1	Gamma
O2_Gamma	0.01031	O2	Gamma

write_csv(bdf, "group_lasso_weights.csv", col_names = T)

And here is the visualization of the corresponding weights as a topomap:

Topomap of Beta weights